A minimal model for chaotic shear banding in shear-thickening fluids 
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We present a minimal model for spatiotemporal oscillation and rheochaos in shear-thickening 
complex fluids at zero Reynolds number. In the model, a tendency towards inhomogeneous flows in 
the form of shear bands combines with a slow structural dynamics, modelled by delayed stress relax- 
ation. Using Fourier-space numerics, we study the nonequilibrium 'phase diagram' of the fluid as a 
function of a steady mean (spatially averaged) stress, and of the relaxation time for structural relax- 
ation. We find several distinct regions of periodic behavior (oscillating bands, travelling bands, and 
more complex oscillations) and also regions of spatiotemporal rheochaos. A low-dimensional trun- 
cation of the model retains the important physical features of the full model (including rheochaos) 
despite the suppression of sharply defined interfaces between shear bands. Our model maps onto the 
FitzHugh-Nagumo model for neural network dynamics, with an unusual form of long-range coupling. 

PACS numbers: 82.70.-y, 05.45.-a, 83.60.Wc 



I. INTRODUCTION 



Complex fluids have long been known to show strong 
coupling between structure and flow. This leads to vis- 
coelasticity in both the linear and the nonlinear response. 
The latter can include both shear thinning and shear 
thickening, with upward and downward curvature respec- 
tively in the steady state flow curve (7(7) of shear stress 
against strain rate. Polymers are usually shear thinning, 
whereas viscoelastic micellar systems can be either thin- 
ning or thickening, as can dense colloidal suspensions 
0, In extreme cases of shear thinning, where (7(7) 
becomes nonmonotonic, steady flow is mechanically un- 
stable on the decreasing portion of the curve. Stability 
can sometimes be restored by shear-banding, in which 
bands of material with unequal strain rate 7 but equal 
stress a coexist, with interface normals in the shear gra- 
dient direction 0,3. For extreme shear thickening, the 
same applies interchanging cr and 7, with band interface 
normals now in the vorticity direction . 

It has recently become clear, however, that in some 
complex fluids (as outlined below) parameter regimes 
exist where the constitutive response to steady driving 
is intrinsically unsteady. This entails a different (or at 
least stronger) dynamical instability to the one present 
in shear-banding. These instabilities are not transient 
phenomena: they rather indicate the presence of com- 
plicated dynamical states within the flow diagrams of 
complex fluids. In principle there can be many sources 
of dynamical instability: under appropriate conditions, 
complex fluids — like simple fluids — are subject to classi- 
cal inertial (Taylor-Couette, turbulence) or thermocon- 
vective (Rayleigh-Benard) instabilities 0, ■ More spe- 
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cific to complex fluids are stick-slip phenomena (like the 
'spurt' or 'melt fracture' effect [ll) and various elastic 
instabilities (so-called 'elastic turbulence' Q)- 

Here we address a distinct class of instabilities in which 
not only the mechanical response, but also the internal 
structural parameters of the fluid, vary in time. Experi- 
mentally, such structural instabilities arise in a variety of 
systems. Typical observations fall into two broad types. 

The first type of unstable temporal behavior com- 
prises sustained, periodic oscillations of the shear rate 
at constant imposed shear stress, or vice-versa: this has 
been observed for surfactant solutions either in worm- 
like micellar plmses [13, IH El IH El or lamellar 
phases 0,0,^3,0], as well as in polymer solutions 
and concentrated colloids ■ The second type of un- 
steady behavior is more complex, with erratic temporal 
responses of either the shear rate or shear stress. Such 
irregular signals have been observed in worm-like mi- 
celles [2ll[il[2ll2l, lamellar (onion) phases [ll[li[23 
and concentrated colloids [2^. In many of these sys- 
tems, such erratic responses occur for parameter values 
enclosed within those where oscillations arise. 

There are strong indications [3, [^ [111 that these 
erratic responses result from a deterministic chaotic dy- 
namics. A remarkable aspect of chaos in these flows is 
that the Reynolds number is virtually zero: the inertial 
term vS/v in the Navier-Stokes equation is negligible. 
The instabilities thus stem from constitutive nonlinearity 
in the rheology of the fluid, unlike the convective insta- 
bility that gives rise to turbulence in Newtonian flows. 
The term 'rheological chaos', or simply 'rheochaos' has 
been coined for such behavior [27l |. 

The involvement of the fluid microstructure in these 
instabilities was established by monitoring structural ob- 
servables and showing that these evolve in concert with 
the time-dependent rheological signal. The methods 
used include birefringence imaging [ill ITp . light scatter- 
ing [HEIl and spatially resolved NMRpi]. Many (but 
not quite all) instances of structural instability occur in 
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shear-stress or shear-rate ranges close to non-equilibrium 
transitions between distinct phases or textures in the 
fluid, for example the transition from isotropic to flow- 
aligned- nematic structures in worm- like micelles '24], or 
the disordered-to-layered packing transition in multil- 
amellar onions [iMIlaillSlS^- (There may also be under- 
lying transitions from a flowing to a jammed state in col- 
loids ,26] : and from isotropic to string- like structures in 
polymer solutions Structural instabilities arguably 

arise when the fluid under flow hesitates between possible 
alternative structures near such transitions. 

A crucial question is whether structural instabilities 
are purely temporal or spatiotemporal in character. Do 
all points in the fluid follow the same time evolution, 
or do different parts of the fluid have different mechan- 
ical and structural states at the same instant of time? 
Early experiments on worm-like micelles IT^ and 
polymers [19| pointed towards heterogeneity: optical ob- 
servations showed alternating turbid and clear bands. 
More recent advances allow spatially and temporally re- 
solved measurements of velocity profiles within a rheome- 
ter. Such experiments, on multilamellar onions 17, 2^ 
and worm- like micelles j24j , unambiguously demonstrate 
both the presence and the nature of the heterogeneities: 
they are fluctuating shear bands. Theoretical models that 
address only temporal instability in a spatially uniform 
system, such as that of Ref. are therefore of limited 
relevance to the experimental situation. 

In this paper we extend the model of Ref. 27' to allow 
for spatial heterogeneity, exploring the resulting scenario 
of why shear bands destabilise in such fluids and how 
this produces oscillatory and chaotic flows. (Note that 
we work in a parameter regime where the flow curve is 
monotonic, so that shear bands cannot exist as a time- 
independent steady state.) 

We thereby create a minimal model of spatiotempo- 
ral instabilities in shear-thickening fluids. The model 
includes an intrinsic short-time tendency to form shear- 
bands coupled to a slow relaxational component dynam- 
ics for the fluid microstructure (modelled as a retarded 
stress response term). In the purely temporal model 
of Ref. 1271 oscillations but not chaos were found (un- 
less a physically unconvincing 'double memory' term was 
used). By allowing for full spatiotemporal dynamics, we 
show that the interplay between the above two factors 
gives rise to several distinct periodic regimes (including 
oscillating shear bands and travelling bands) and also 
regimes of spatiotemporal rheochaos. Preliminary ac- 
counts of our work were given in Refs. lsollsil 

Our motivation for studying the case of a fluid that 
shows shear-thickening (in itself a wi dely observed but 
poorly understood phenomenon |2^, [23) is that sev- 
eral of the above-cited experiments concern such fluids 
(Refs. ES mill El I23 for wormlike micelles; [l^ for 
polymer solutions; and |2(ll2^ for colloidal suspensions). 

The alternative case of shear-thinning fluids has been 
recently addressed by Fielding and Olmsted 's^, and in 
less detail also by ourselves ISlj- The authors of Ref. 32 



have demonstrated the presence of a rich dynamics, in- 
cluding rheochaos, in their model. Also closely related 
is the work by Chakrabarti et al. [s^l on nematic liq- 
uid crystals which also shows regimes of spatiotemporal 
chaos. In what follows, we shall note similarities and dif- 
ferences between our own work and these other studies. 

The rest of this paper is organized as follows. In sec- 
tion^ we define the model, discuss its physical assump- 
tions, and relate it to other models (both in rheochaos 
and in other fields of study). In section UTTl we give a 
qualitative description of how the model works. Sec- 
tion IIVI explains how the model is solved numerically. 
In section we present our results in the form of a 
nonequilibrium 'phase diagram' and comment on the var- 
ious flow regimes encountered. Section IVII studies a low- 
dimensional version of the model which offers further in- 
sights into the nature of rheochaos. In Section rVIII we 
summarize our findings and their generic implications for 
the physics of rheological instabilities in complex fluids. 



II. THE MODEL 

Our model has two main physical ingredients, as fol- 
lows: 

(i) the fluid has a tendency to form shear bands; 

(ii) there are slow structural modes in the fluid whose 
delayed relaxation modifies the evolution of the stress. 

The latter ingredient is central to our description. 
When some structural mode is disturbed by the flow, 
it will relax on a timescale that is distinct from the stress 
relaxation time in the system. Below, for simplicity, we 
will consider a single timescale rs for structural relax- 
ation, with stress relaxation assimilated into the usual 
Maxwell time tm- We focus on the case where structural 
modes relax at least as slowly as the stress itself would 
do at fixed structure (rs > tm). 

Because of the disordered energy landsca pe and 
metastability intrinsic to many complex fluids '3^ l35l 
l3^ . such slow dynamics is commonplace. Candidates for 
slow structural modes include the mean length or the lo- 
cal gel fraction in worm-like micelles; local composition 
variables (e.g., colloidal volume fraction); and 'fluidity' 
parameters reflecting, e.g., a local bonding state |3lll37l|. 
An involvement of slowly evolving structure has been 
shown in many experimental cases 0, IT3 . 

ElESIilIillii. 



A. Model equations 

We consider the situation of a fluid under pure shear, 
and we assume, as usual 0, that the shear stress a de- 
couples from other stress components, and depends only 
on the rate of shear strain 7. We will restrict our study 
to one-dimensional spatial heterogeneity in the fluid, in 
the vorticity direction (perpendicular both to the velocity 
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FIG. 1: Vorticity shear bands in a Couette cell 0,0. These 
occur in shear-thickening materials with an S-shaped flow 
curve (a vs. 7). A given mechanical torque on the Couette im- 
poses a mean (spatially averaged) stress (cr) . Any (cr) within 
the unstable portion of the flow curve is unstable toward shear 
bands (with different local microstructures) along the vortic- 
ity direction z, depicted here as a clear and a turbid band. 
These coexist at a common shear rate 7c (whose value is flxed 
by gradient terms Q) but at different stresses a\ and (T2. The 
amount of each band is such that the weighted mean of a\ and 
(T2 matches the externally imposed mean (cr). 

and the velocity gradient); in cylindrical Couette geom- 
etry, which we will always refer to in the following, this 
corresponds to the axial coordinate denoted z. (For 3D- 
related effects in situations of shear-banding, see 38].) 

This choice for the heterogeneity direction is motivated 
by the classical geometry of steady shear bands in thick- 
ening materials where, as shown in Fig.^ bands are 
stacked one onto another in the vorticity direction (vor- 
ticity shear bands). Note that, under these assumptions, 
the shear flow is homogeneous within each slice of height 
z, as required by the low- Reynolds limit. (We neglect any 
small variation in the velocity gradient direction caused 
by the curvature of the cell. Without fluid inertia, the 
stress cannot vary in this direction.) 

In all of the following, we shall work under conditions 
of imposed torque, i.e., under an imposed value of the 
mean (spatial average) of the stress (a) : this is the usual 
situation for vorticity shear bands, shown in Fig. [We 
also studied the behaviour of our model under condi- 
tions of imposed shear rate (results not shown) but no 
spatiotemporal effect was observed, see note 

In our model, the shear stress cr(z,t) at time t evolves 
according to the following equation of motion: 

ct(z, t) = 7(<) - i?(cr) - A I M(t-t') a{z, t') dt' + kV^ct 

J — 00 

(1) 

where 

R{a) ^ aa - ba'^ + ca^ and 7W(u) — 6""/^^ (2) 

Eqs. ai'o non-dimensional: the transient elastic 

modulus is taken as the unit for stress, and if, the axial 
extent of the Couette cell (vertical in Fig. ^ , the unit of 
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FIG. 2: Plots of R{a) and the inverse function a{R). The 
S-shape of the latter encodes the tendency of the fluid to form 
vorticity shear bands. 

length. Several choices are possible for the time unit; the 
simplest one is the fluid's Maxwell time tm- However, 
as we will always be dealing with times longer than tm, 
we choose instead 100 tm as a time unit. [This amounts 
to setting tm = 0.01 and, via Eq. (PI, a = 100 through- 
out what follows; we do this.] The reader is referred 
to Appendix 1X1 for details on the non-dimensionalisation 
procedure. 

Note that & ~ da/dt in Eq. © is a local time deriva- 
tive. We also emphasise that the shear rate 7 in Eq. ^ is 
uniform: the moving wall of the rotor imposes the same 
velocity for all heights z, and 7(2, t) — j{t) only. 

The terms in Eqs. and |2Jl have the following sig- 
nificance. i?(cr) describes instantaneous, nonlinear stress 
relaxation. As defined by Eq. (0), is a third-order poly- 
nomial, with the positive constants a, b and c chosen so 
that R{(j) > for cr > 0, and so that the inverse function, 
cr{R), is an S-shape (see Fig. We choose a — 100, 
b — 20, c = 1.02 in this paper; setting c = 1 would 
give a physically inappropriate zero of i?(cr) at cr = 10. 
The i?(cr) term encodes ingredient ft) of the model, i.e., 
the tendency of the fluid to form vorticity shear-bands. 
This tendency is instantaneous, but frustrated by struc- 
tural relaxation effects — see Fig. 13 Note that, with our 
choice of units, a controls linear stress relaxation in the 
fluid, and hence determines the Maxwell time tm' 

Tm = 1/a (3) 

Setting & = c = A = K = Oin Eq. indeed recovers the 
classical Maxwell model for linear viscoelasticity. 

The integral term in Eq. © corresponds to our second 
physical ingredient, and represents retarded stress relax- 
ation due to slow structural reorganisation in the fluid; A 
is a positive constant governing the strength of this term. 
For the sake of simplicity, retarded relaxation is chosen 
linear in past stresses; but note that, as discussed in 
Ref. I27I a similar form can be obtained by introducing an 
explicitly structure-dependent i?(cr) and then linearizing 
this structure-dependence. The memory kernel A4 is in 
principle any decaying function p^ . However, in Eq. I©, 
we choose it mono-exponential with characteristic time 
Ts. This choice permits a much simpler, fully differential 
representation that we exploit in our numerics below. As 
stated previously, we take structural relaxations to be 
slow compared to the intrinsic time scale of stress relax- 



4 



ation; we study mainly 4 < ts/tm < lO**. (This contrasts 
somewhat with the model of Ref . Is^ where ts/tm — 1-) 

Finally, the last term of Eq. ^ assigns to stress a dif- 
fusivity K. By analogy with the classical Cahn-Hilliard 
model for phase separation |33j, this can be viewed as 
representing the cost of maintaining interfaces in inho- 
mogeneous states. Such minimal nonlocality in the con- 
stitutive model is known to be physically crucial. For 
example, in steady-state shear bands it drives a selection 
mechanism, among all possibilities on a given flow curve 
like that of Fig.^ for the coexistence strain rate 7c Q. 

A superficial comparison of Eq. Q with Ref. [23 might 
suggest that our work is new only in the addition of this 
diffusive term. But of course we simultaneously upgrade 
the dimensionality of the problem from a{t) in that work 
to (t(z, t) here; this allows vastly richer behavior, in better 
correspondence with the experimental facts, to emerge. 

We now recast Eqs. Q"® a purely differential 
form Q^. We define a new variable 

m{z,t)^ f M{t~-t')a{z,t')dt' (4) 

J — QO 

which we call the 'memory', and which contains the struc- 
tural part of the stress relaxation. We can then rewrite 
Eq. Q as an exactly equivalent differential system 

(T = 7 - R{a) - Am + kVV (5) 



Equation (jSJ states that the memory is always relaxing 
towards the current and local value of the stress cr{z,t) 
with a (slow) rate t^^. 

In its differential form, the model can be interpreted as 
follows. Equation © is a nonlinear rheological equation 
controlling the temporal evolution of the stress, where 
one of the participating variables is structural by nature 
(the memory m); this structural variable is subject to 
a distinct dynamics, governed by a 'structural equation', 
Eq. © , with a different relaxation time rg . The coupling 
of the two dynamics, mechanical and structural, is the 
source of the instability. 



B. Relation to other models 

Viewed in this way, our model belongs to a larger class 
of models that have been recently proposed to describe 
the dynamics of complex fluids. 

The models of Refs. [3^ and [U for shear-thinning 
solutions of worm-like micelles work on a similar scheme: 
there, the 'structural variable' is the mean chain length of 
the micelles ^| , or simply the Maxwell time of the fluid 
itself 'sij. In both cases, its evolution is governed by a 
differential 'structural equation' akin to Eq. Belong- 
ing to the same family, albeit coming from a somewhat 
different perspective, Derec et al. proposed a model for 



the fluidisation transition in pastes |37|, where a classi- 
cal rheological equation is supplemented by a structural 
equation describing a 'fluidity' parameter. 

From a more formal point of view, all these models 
are related to the prey-predator and reaction-diffusion 
models commonly found in nonlinear physics and math- 
ematical biology (III our case, a and m are the 
two competing species). This type of non- linear models 
are known to yield complex spatiotemporal behavior j 
thus it should not come as a surprise that the related 
models used in the field of complex fluids do also. 

It is interesting to note that our Eqs. 0"® in 
fact map exactly onto a well-known model of nonlinear 
physics, the FitzHugh-Nagumo model ^Mi,^M,M- 
The simplest version of the latter model was developed to 
describe electrical activity in a single neuron, where the 
axon membrane potential is a fast variable (also called 
'excitation variable', analogous to a in our model), and 
is coupled to the dynamics of a slow variable (or 'recov- 
ery variable', analogous to m), related to the activity of 
sodium ion channels. This system is known to produce 
van-der-Pol-like (purely temporal) oscillations, like those 
described for complex fluids in the earlier, spatially ho- 
mogenised, version of our model 0|. 

In recent years, spatially inhomogeneous extensions of 
the FitzHugh-Nagumo have been used to study the col- 
lective properties of interacting networks of neurons, and 
other related assemblies of coupled nonlinear oscillators, 
focussing on the competing effect of local and global cou- 
plings miiHilliil. The model of Eqs. ©-® is of just 
this type: the local coupling is supplied by the diffusion 
term kV^ct, and the global coupling by the constraint 
that the spatial mean stress (a) is externally fixed. 

We remark however that: (a) in our model, the mean 
value of the stress is fixed, while in neural networks where 
a coupling to the mean membrane potential is imple- 
mented jlH . W\ li^ . its value fluctuates in response to 
the collective dynamics in the model; (b) possibly for 
this reason, several of the spatiotemporal patterns re- 
ported below (in particular oscillatory shear bands) have 
not to our knowledge been reported in the literature on 
FitzHugh-Nagumo networks; (c) an active research topic 
in neural networks is the effect of noise on global behav- 
ior, including the possibility to trigger chaotic dynam- 
ics . Such a role for noise (either thermal, or mechan- 
ical) is neglected in most models of rheochaos but would 
constitute an interesting avenue for future research. 



III. QUALITATIVE FEATURES 

We now explain some qualitative features of our model, 
focussing on the origin of the dynamical instability, and 
on a physically important scaling property. We also 
briefly discuss how the values of our model parameters 
may be estimated from experiment. 
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FIG. 3: Steady-state flow curve (thick line), and underlying 
short-term flow curves (thin lines, with, from left to right, 
m — 0,1,2, ... , 20). The stress range between the dotted lines 
corresponds to the unstable region, here deduced for rs = 100 
from Eq. 0. Numerical parameters: A = 40, a = l/T-M = 
100, 6 = 20, c = 1.02, K = 0.01, H = 1. 



A. Flow curves 

We start by computing the flow curves for the fluid, 
i.e., the relation between a and 7 in steady-state flow. 

From Eqs. the flow curve for steady-state, ho- 

mogeneous flows is given by setting & = m = V^cr = 0. 
Equation Q then yields m — a: at each point on the 
curve, the memory has relaxed to the steady-state stress. 
This can be substituted in Eq. ((SJ, to provide the equa- 
tion for the steady-state flow curve |27|: 



7 = R{a) + Act 



(7) 



When R'{a) + X > for all values of a (a condition which 
we shall always take to hold below) the flow curve is 
monotonically increasing, as in the thick curve in Fig. |31 
In contrast with rheological constitutive models that do 
not show structural coupling, this monotonicity does not 
ensure mechanical stability of homogeneous flows 27] . 

The qualitative reason for instability of the monotonic 
flow curve is that the steady states associated to this 
curve only arise when the memory m has relaxed, i.e., 
for time scales beyond Tg. For times much shorter than 
rs, the fluid instead behaves as if the memory m was 
frozen. Thus, at short time scales, the fluid lives on one 
of a set of 'instantaneous flow curves' for which stress has 
relaxed, but not structure. These instantaneous curves 
have (7 = at flxed to, giving 



7 = Ric) + Am (short term). 



(8) 



where m is a parameter. Such curves for to = 0, . . . , 20 
are plotted in Fig. O and are not monotonic. As the 
memory slowly relaxes, the fluid drifts from one such 
curve to another; Eq. (jT)) can be reconstructed by pick- 
ing, for each value of the stress cr, the corresponding point 
on the particular instantaneous curve that has m = a. 
But for our choice of non-monotonic R, the fluid has an 



instantaneous tendency to form shear-banded flow. This 
impedes the establishment of a steady state and insta- 
bility can arise even when the steady-state flow curve is 
monotonic 123 • 



B. Origins of the dynamical instability 

As just discussed, the existence of decreasing portions 
of the short-term flow curve, in contradiction with the 
steady-state flow curve, is a source of instability. In the 
homogeneous version of the model (no z-dependence) , 
this causes temporal instability in the form of van-der- 
Pol-type oscillations |23| . Essentially, the mechanism for 
such oscillations is as follows: starting for instance from 
a situation of high stress as compared to the memory 
(cr > to), the memory m will start to grow as dictated 
by Eq. ©; this growth of to in turn brings a decrease 
in the value of a through Eq. (O; adapting to this, to 
then decreases, which increases cr, thereby recovering the 
initial high-stress situation and starting the cycle anew. 

When spatial dependence is added, this basic tempo- 
ral oscillation becomes compounded with shear banding. 
It is easily seen how this can give complex spatiotempo- 
ral dynamics: looking at Fig. 13 if one imposes a mean 
stress (cr) chosen within the unstable region, the system 
decomposes into several shear bands with unequal local 
stress, as depicted on Fig. ^ Unlike the classical situa- 
tion of Fig. here the van-der-Pol-type temporal oscil- 
lation rules out steady states for the shear bands: a local 
unstable dynamics engages. Simultaneously, the bands 
have to match, all together, the constraint on (cr). This 
creates long-range couplings between the bands (when 
one oscillates, another has to compensate), considerably 
complicating the spatiotemporal behavior. 



C. Linear stability analysis 

We have performed a standard linear instability calcu- 
lation on the model 0-®, with the following results: 

(i) for an externally imposed value of the mean 
stress (cr), instability for stress evolution arises only if 



R'{{^)) + — < 0. 



(9) 



(a) when the above instability condition is obeyed, 
only stress modes with wavevectors q in the range 



< g < 



R'{{a)) + l/rs 



(10) 



are unstable. Also, modes with smaller wavevectors have 
a larger growth rate: when the uniform q — Q mode is 
constrained, the lowest nonzero g-mode grows fastest. 

The instability condition (i) is the same as for the 
purely temporal version of the model |23|; it states that 
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instability can arise only in regions of decreasing R (or 
equivalently, within the decreasing portions of the short- 
term flow curves in Fig. , and that it is facilitated as 
Ts becomes larger or impeded as it becomes shorter. 

The second result (ii) is a natural consequence of stress 
diffusion with diffusivity k; this kills off fluctuations with 
too large q (small wavelengths) and prevents such modes 
from becoming unstable. 

D. Scaling property of the model 

Let us assume that a stress response a{z, t) and shear 
rate j{t) are obtained as the solutions of the model of 
Eqs. 151)-®, under conditions of imposed mean stress 
(a) and with a given set of parameters a = Tj^^, h, c, 
A, TS, K. Then, choosing any scaling factor a, the scaled 
responses aa{z,t) and a'^(t) would also be obtained as 
the solutions of the model if one applied a mean stress 
a{a) and used a scaled set of parameters a = Tj^^, 6/a, 
c/a^. A, Ts, K. 

The proof of this scaling property is elementary and 
left to the reader. It means that, although all numeri- 
cal results given below are found with one specific set of 
model parameters (A = 40, a = I/tm = 100, b — 20, 
c = 1.02, K — 0.01, H = 1), they in fact describe the 
behavior of a one-parameter manifold of parameters gov- 
erned by the scale transformation parameter a. 

Secondly, for our chosen parameter set, the observed 
stresses in both oscillatory and chaotic regimes are very 
large (of order 10 in dimensionless units, i.e., ten times 
the elastic modulus of the fluid). While such large am- 
plitudes of the elastic strain may be relevant for certain 
types of fluids, no direct importance should be attached 
to this: one can scale down numerical values of the stress 
to around, say, unity, keeping exactly the same oscilla- 
tory or chaotic dynamics. Note however that this will 
not change the ratio of stress fluctuation to mean stress; 
this is always large for our parameter settings. (The pa- 
rameter space is vast, and we have not explored values 
for which this ratio might be reduced.) 

E. Experimental estimates for model parameters 

Our model involves a set of phenomenological model 
parameters whose numerical values are a priori unknown. 
One way around this difhculty is try to extract esti- 
mates for these values using appropriate experiments: it 
is thus useful, at this point, to sketch some of the possible 
strategies — keeping in mind, however, that our minimal 
model is intended essentially for a qualitative exploration 
of the physics at work in complex fluids not for an effec- 
tive fitting tool for specific nonlinear materials. 

(It is also recalled that numerical values used in the ar- 
ticle are non-dimensional; see Appendix^for the proce- 
dure to convert dimensional quantities, as obtained from 
experiments, to non-dimensional ones.) 



A first experimental strategy is as follows: starting 
from a state of complete rest (where the structural mem- 
ory has fully relaxed, m — 0), and assuming that the 
structural time ts is slow enough, one may be able to 
measure the instantaneous flow curve at zero memory 
[Eq. (jS))] using a moderately rapid ramp up of the ap- 
plied stress (e.g., sample the whole curve within a few 
minutes). Of course, the obtained curve would be unsta- 
ble (fig. Oil, showing a conventional 'discontinuous shear- 
thickening' scenario with a shear-rate plateau and classi- 
cal (steady-state) shear bands; but the information col- 
lected on the plateau (critical shear rate, lower and upper 
values of the stress at the jump, . . . ) would probably be 
sufficient to estimate the short-term parameters a, 6, c, 
and thereby determine the fimction R{(t) in the model. 

Following this determination of R, the value of the cou- 
pling parameter A could then be found through a simple 
steady-state experiment where a constant stress CTss is 
imposed for a long period a time (in a stable region of 
the phase diagram): as stated by Eq. 10, the shear rate 
then corresponds to the (known) instantaneous contribu- 
tion i?((Tss) plus the delayed contribution Actss- Thus, a 
measure of the shear rate will determine the value of A, 
the only unknown in the equation. 

Next, estimating the value of the structural relaxation 
time Ts could be done through a relaxation experiment: 
from a situation of steady state under a given applied 
stress, the value of the stress is suddenly changed to a new 
constant value (or even cancelled) ; if ts is long compared 
to other times in the system, it should be controlling the 
global relaxation of the response to the new steady state. 
Another possibility which could be explored could be to 
look at the linear viscosity after a period of strong shear 
and/or after a temperature jump j29j. 

The last parameter that needs to be determined is the 
stress diffusivity k, which in classical, steady-state shear 
banding is related to the width of the interfaces between 
bands. We expect k to be extremely small ^3 (and 
in any case, much below the artificially large value of 
10~^ used throughout this work for computing reasons), 
although we are not aware of any experiments allowing 
to measure it. This should not be too much of an issue, 
as our results show that (below a certain threshold) the 
behaviour of the model does not qualitatively depend on 
the actual value of k. 



IV. NUMERICAL METHODS 

In this section, we detail the numerics techniques used 
to obtain the results presented in the rest of the article. 



A. Spectral scheme; boundary conditions 

Our numerical scheme expands the stress in Fourier 
modes, then is solving for the evolution of these — rather 
than directly solving the model equation on a grid of 
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points in direct space |33,|33|- Such a spectral scheme |50j 
had two main advantages, as follows: (a) in Fourier 
space, working under conditions of fixed mean stress sim- 
ply corresponds to fixing the value of the uniform, zeroth 
mode in the Fourier expansion of the stress (this is much 
more difhcult to implement in direct-space schemes); 
(h ) Not only can arbitrarily accurate numerical results be 
obtained by keeping enough Fourier modes in the scheme 
('high-order truncation'), but also, by keeping only a 
minimal number of modes ('low-order truncation'), one 
can obtain a reduced description whose analysis is much 
simplified. Such a truncation is more likely to capture the 
physics of the problem than a reduced real-space scheme 
with a spatial grid of only a few points. We shall present 
both high-order and low-order numerics in what follows. 

We must also choose appropriate boundary conditions 
on the system. We demand zero stress flux at both ends 
of the Couette cell: 



Vcr|z=o = V(T|z=ff = 



(11) 



To motivate this, note that stress flux arises through the 
diffusion term in Eq. (jSJ, which is generally ascribed to 
material displacement of stress-carrying elements |5l| . 
Our boundary conditions then reflect the fact that no 
fluid element leaves through the top nor bottom of the 
cell. 

Now recall that, for reasons related to the physics of 
vorticity shear banding (Fig. we have the constraint: 



1 

n 



H 



a{z, t) dz = const. 



(12) 



We now decompose the stress field onto a Fourier basis 
of spatial modes compatible with the above boundary 
conditions (cosines only): 



a{z,t) ~ ^ ak{t)cos{kTrz/H) 



(13) 



k=0 



and note that ao{t) = (a) (given that iJ = 1 in our 
units). Similarly for the memory term: 



to(z, t) 



E 

k=0 



mk{t) cos{kTr z/H) 



(14) 



The series in Eqs. (|13f) and (|14|l only contain the first N 
modes ak{t) and mk{t) of an infinite Fourier series: the 
level N of this truncation (known as a Galerkin trunca- 
tion [^) determine the overall accuracy of the scheme. 
(For a given accuracy, the number of modes that must be 
kept is often much lower than the number of grid points 
in any real-space discretisation 50J.) 

We now project Eqs. (|13|I -H14 |) onto the Fourier ba- 
sis to obtain a set of evolution equation for each of the 
mode amplitudes ak{t) and mk{t). Thanks to the sim- 
ple functional form of the model, it is possible to obtain 



analytical expressions for the mode equations, which are 
given in full in Appendix^ These mode equations are 
of the form (1 < fc < - 1) 

(Tfc — (coupling terms) — Xnit — Kql<Jk (15) 
nik - (7k 



rrik 



TS 



(16) 



where qk = kir/H is the wavevector, and the coupling 
terms, stemming from the non-linearity of the instanta- 
neous term R{a), link the evolution of each ctj, mode to 
all others. As expected, stress diffusivity damps higher 
modes via a linear Kq^ term. 

Finally, the equations for the fc = uniform modes 
cro(t) = (cr) and mo(t) are 

1 

ao = O^iit)-- J R{a{z,t))dz-Xmoit) (17) 
((t) - m,o{t) 



mo 



(18) 



Integrating Eq. l(TS|l gives [with m{t = 0) = 0, see below] 



mo(t)-(a)(l-e-*/^^) 



(19) 



Plugging this expression into Eq. p7|l . we obtain an im- 
portant expression which provides us with the instanta- 
neous value of the shear rate 7, once the value for the 
stress has been calculated from the Fourier modes: 



1 



H 



iit) = j^ RH^. t)) dz + A(a) (1 - e-*/^^ ) (20) 

The spatial integral in this expression can in fact be car- 
ried out analytically, not numerically (see Appendix IB|) . 

The final outcome of the Fourier-Galerkin truncation 
scheme is that the partial differential equations (0-® 
in the original formulation of the model have become a 
dynamical system of finite order, containing 2A^ ordinary 
differential equations for the modes ak{t) and mk{t) as 
specified by Eqs. (|15|l - Hlt)|l . plus the above equations for 
the uniform modes. The numerical integration of the 
ordinary differential equations (|15|I - H16|) was performed 
with the help of a commercial solver package |52 | , using 
an adaptive time step. Given the separation of scales 
between the typical times tm and Tg, the adaptive step 
was needed, to maintain acceptable computing times. 

Initial conditions aXt — Q were chosen as follows: mem- 
ory was set to zero (i.e., mfe(O) = for all k), and initial 
values for each stress mode (Tfc(O) were picked at random 
between and 10^^. The noise amplitude was found not 
to be essential to the qualitative results obtained. 



B. High-order and low-order truncation 

The behavior of our fluid model was explored using two 
types of Fourier-Galerkin schemes: one where the number 
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of modes kept in Eqs. is high {N = 40), and 

one where it is minimal {N = 3). 

For the high-order truncation (Section 0, we found 
that keeping iV = 40 modes was numericahy accurate, 
with reasonable computing times. One criterion in this 
choice was that numerics should be able to fully resolve 
interfaces between shear bands (the sharpness of inter- 
faces is essentially controlled by the stress diffusivity k; 
we set it to 10^^ unless otherwise stated). A second 
criterion is that TV should be high enough that results 
become TV-independent: this was indeed the case, with 
truncations from iV = 25 up all giving similar outputs. 
This applies at the level of phase diagrams, etc., but not 
to individual trajectories which (at least in chaotic re- 
gions of the phase diagram) can depend on every detail 
of the numerics. 

We also checked the validity of our high-resolution re- 
sults, for a given set of parameters, against variations in 
either initial conditions for the stress modes, or in the 
actual sequence of time steps followed by the adaptative 
iterator during numerical integration. We found a lim- 
ited dependence on these factors from one run to another. 
To take this into account, the 'phase diagram' of Fig. 01 
where the general behavior of the model is summarized, 
was established on the basis of several independent runs 
for each point marked in the figure. Thus the general 
features of the phase diagram, as well as all conclusions 
drawn on the model's global behavior, are reliable. 

In the low-order truncation fSection lVip . on the other 
hand, our aim was to study the behavior of the model in 
its simplest possible spectral representation. In particu- 
lar, we are interested in whether low-dimensional chaos 
was present and what route led to it. Since three de- 
grees of freedom are required to allow for chaos, the low- 
est compatible truncation in our model corresponds to 
A'^ = 3: excluding uniform modes (Tq = (a) and mo, which 
are dynamically inactive, this leaves a four-dimensional 
system involving cti, a2, rrii and TO2. In the following 
sections, we present the results of the model, first in the 
high-order truncation, then in the low-order. 

V. HIGH-ORDER RESULTS 

In this section, we present the results that were ob- 
tained by solving the model of Eqs. H15|) - (|16|l in the high- 
resolution truncation {N = 40). 

A. Phase diagram 

To get a global picture of the model's behavior, we 
made a systematic exploration of its spatio-temporal dy- 
namics by varying the two physically most important 
parameters: the ratio between the structural timescale 
and the mechanical timescale ts/tm, and the spatially 
averaged stress (a) (fixed by the mechanical torque on 
the Couette). On varying ts/tm, we choose to keep tm 



fixed (tm = 1/a = 0.01) and vary ts only. All other pa- 
rameters are held at their values of Fig. O i.e., b = 20, 
c = 1.02, A = 40, K = 0.01, H = 1. 

Despite the high-dimensional dynamical system under 
consideration (2 A = 80), the obtained 'phase diagram' 
for the model (shown in Fig. ^ displays a simple overall 
structure where three main dynamical regimes emerge: 
periodic response with (more or less complex) oscillat- 
ing shear bands at extremely long ts; periodic response 
with travelling bands for long Tg; and finally, chaotic re- 
sponse at shorter rs and off-centered values of (ct). The 
dotted lines separating regions on the phase diagram 
are guides to the eye, representing crossovers not sharp 
transitions between different types of behavior. The 'C- 
regions' marked in the plot were defined so as to en- 
close all observed chaotic points; these regions do, how- 
ever, contain internal structure with periodic and chaotic 
pockets, whose exact boundaries can depend on initial 
conditions and other details. 

We next discuss in some detail the three main regimes 
encountered on the phase diagram, emphasising an in- 
tuitive understanding of the physics involved. However, 
as shall be seen, there is significant variety in behaviour 
even within the same regime. 



1. Periodically oscillating shear bands 

We first discuss the regime encountered when struc- 
tural evolvution is much slower than stress relaxation 
(typically ts/tm > lO'^) marked 'O' in Fig. ^ Here peri- 
odic oscillations of the shear rate and stress are observed, 
while spatially, the flow presents oscillating shear bands. 
Depending on the imposed stress (a), the waveforms of 
these band oscillations (which induce oscillations in both 
the local stress a{z,t) and the mean strain rate 7(i)) 
range from simple to very complex. Moving along any 
horizontal line within the region O (changing {a) at con- 
stant ts/tm), one observes the same succession of behav- 
iors; we present these for one typical line, ts/tm = 10"*. 

a. 'Flip-flopping' shear bands. Near the middle of 
this line, e.g., {a) = 7.0, we find two 'flip-flopping' bands 
(see Fig. |SJa) : at any time, the cell is equally divided be- 
tween a high-stress and a low-stress band, but the iden- 
tity of the bands reverses periodically. This results in the 
'checkerboard' spatiotemporal pattern shown in Figure[Sl- 
a. Accordingly, the stress a{z,t), measured for a given 
height within the cell (e.g., z — 2/3), displays periodic 
oscillations with a waveform close to a square wave (and 
abrupt changes between the low- and hi gh-s tress states, 
as expected for 'relaxation oscillations' |53). The flip- 
flop period TfHp is of order the structural time rg. The 
shear rate 7(i) also presents a periodic evolution, albeit 
with a more complicated waveform. (This extra com- 
plexity is generic in our shear rate time series, but seems 
to be a model-dependent feature rather than connected 
to any deep physics.) 

The mechanism underlying the dynamics of flip- 
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FIG. 4: Phase diagram of the model when rs and (a) are varied: (a) chaotic point, (•) periodic point. Three main regimes 
are observed: (O) oscillating shear bands, (T) travelling shear bands, (C) chaotic regions. The outer dashed line is the linear 
stability limit, R'{{a)) + 1/rs = [Eq. Q]. Numerical parameters as in Fig.^ 



flopping bands is as outlined in Sec. IIII Bl because of 
the short-term instability present in the fluid (Fig. |2Jl, 
shear-bands form. Each of these bands is locally submit- 
ted to a van-der-Pol-type instability which forces them to 
oscillate between states of high and low stress. But these 
local oscillations cannot occur independently: because 
the spatial mean of the stress (a) must be conserved, 
they have to be synchronous — when one band goes up, 
the other must flip down. Note that shear bands not 
only display different values of the stress, but, in general 
also of the memory m(z,i), albeit much less markedly: 
in our model, this means that both the mechanical and 
structural states differ between bands. 

Finally, a physically important question in this regime 
is: what happens to the interface between bands during 
flip-flops? Looking at Fig. [SJ-a, we see that the interface 
position seems stationary. A fixed position for the inter- 
face would imply that, in the short interval when bands 
are flipping, the interface proflle has to reverse its slope. 
Interestingly, closer scrutiny shows that this is not what 
happens: instead of reversing the slope of a static inter- 
face, the system prefers to replace the 'wrong' interface 
by a new one with the correct slope. As seen in Fig. El 
this occurs by a rapid sweeping motion where the old 
interface quickly travels to one end of the cell and disap- 
pears; simultaneously, the new interface is generated at 
the other end and moves to the middle of the cell, where 
it will rest until the next flip occurs. Thus, long peri- 
ods of immobility for the interface alternate with rapid 
sweeping motions (the latter being too quick to be ob- 
servable on the scale of the plot in Fig.|5l-a). 

Qualitatively, everything happens as though when flips 



occurs in the bands, they trigger a travelling wave remov- 
ing the old interface and bringing the new one in. Thus, 
this regime of flip-flopping bands must be seen as one of 
intermittent travelling waves, separated by long periods 
of latency (of order ts/2). In accord with this view, when 
the latency interval diminishes sufficiently (i.e., rs is de- 
creased), the intermittent waves should become continu- 
ous: this is indeed exactly what happens in the 'T'-region 
of the phase diagram (discussed in Section IV A 2|) . 

b. Oscillating bands with zig-zagging interface. Pur- 
suing our exploration of the model's behavior along the 
horizontal line ts/tm — 10*, we now move slightly off- 
center to the right or the left; for example, (ct) = 7.1. 
The interface between bands now adopts a periodic, zig- 
zagging motion which superposes to the synchronous os- 
cillations of the shear bands (see space-time plot in Fig.jSj- 
b). Accordingly, the time series of the stress becomes 
more complex (square waves are distorted), and the pe- 
riod becomes a multiple (here three) of Tflip. A similar 
increase in complexity is seen in the shear rate. 

The zigzag motion arises because the off-centered value 
of (cr) now enforces unbalanced proportions of the low- 
and high-shear bands in the Couette cell. Because this 
value for (cr) is fixed, these proportions are also fixed: 
thus each time shear bands oscillate and reverse identi- 
ties, the interface between them must move to and fro to 
maintain the required proportions. 

c. Complex periodic oscillations. If we now move 
further towards the wings of the phase diagram, the 
coupling between the global constraint on (cr) (unequal 
bands) and the intrinsic flipping dynamics of the bands 
generates an extremely complex behavior in the fluid. 
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FIG. 5: Typical periodic responses in the oscillating shear 
band regime, chosen along the line ts/tm = lO** in the phase 
diagram, (a) Flip-flopping bands at (a) — 7; (b) Zig-zagging 
interface at (a) = 7.1; (c) Complex periodic motion at (a) = 
9. Each group presents time series of the stress ct at « = 
2/3, the shear rate 7, and a space-time plot of cr(z,t) with t 
vertical, z horizontal (clear shades correspond to high stress, 
dark shades to low stress). Parameters as in Fig. |21 
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FIG. 6: Plots of (j{z,t) and m(z,t) vs. 2 for successive times, 
showing how the 'old' interface is replaced with a 'new' one in 
the regime of flip-flopping bands. In (a) and (b) one sees the 
latency period with stationary interface — note however that 
m slightly evolves from (a) to (b), thus preparing the forth- 
coming flip; (c) As the flip occurs, travelling wave replacing 
the previous interface by a new one; (d) and (e) new latency 
period — note again the slow evolution in m. 



which, remarkably, manages nonetheless to remain peri- 
odic. Figure [S}c shows the response obtained at (<t) ~ 
9.0: the time series for the stress (t(2/3, i), though still 
related to the simple oscillations seen at (a) — 7.0, has 
become very ragged, and the period is now a large mul- 
tiple of Tflip (about six times). 

We conclude this presentation of regime O by not- 
ing that, in experiments, oscillating shear bands similar 
to our theoretical results have indeed been observed in 
shear-thickening fluids [H IH El El . 



2. Travelling shear bands 

We now describe the second regime in the phase dia- 
gram of our model, encountered when tq is long as com- 
pared to Tjvi, but not exceedingly so (10 < ts/tm < 10^; 
region 'T' in Fig. QJ- This regime is characterised by 
a periodic nucleation of shear bands which subsequently 
cross the system at roughly constant velocity. The typ- 
ical situation is shown in Fig. nucleation occurs at a 
boundary, and bands travel across the entire system, one 
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FIG. 7: Periodic time series and space-time plot of the stress 
c^Zjt) in the travelling band regime. Parameters: tq/tm ~ 
90, (cr) = 7, others as in Fig. El 



at a time. As before, the period of the various time se- 
ries is comparable to the structural timescale Tg . Similar 
band motion was reported in Ref. 32, but there the bands 
'ricochet' off the walls of the container. 

Figure |S1 shows examples of other travelling behaviors 
in this regime: The first space-time plot (a) presents a 
case where nucleation occur at a non-boundary point, 
yielding two outgoing waves in a one-dimensional ana- 
logue of the classical 'target patterns' seen in chemical os- 
cillators [s^l ■ Space-time plot (b) presents another case of 
interior nucleation, but where bands alternatively travel 
to one edge or the other. (Details here strongly depend 
on the initial conditions.) Finally, space-time plot (c) il- 
lustrates how the model behavior crosses over smoothly 
from oscillating bands in regime O to travelling bands 
in regime T: the behavior is intermediate between the 
situations depicted in Fig. [S}a and Fig. [7| with bands 
somehow oscillating as they travel. 

Note that the waves observed in regime T are kine- 
matic waves, arising from a staggered phase distribution 
in the local van-der-Pol band oscillations. They are not 
associated with transport of material or stress; hence the 
waves can travel through the system boundaries, despite 
the imposed boundary conditions (Vcr = 0). Accordingly 
the wave velocity is independent of the stress diffusion 
constant k (data not shown). 

Finally, we studied how the bands velocity varied as rs 
was changed. At least for the situation in Fig. [7| where 
there is only one band at a time in the cell, a naive argu- 
ment for the velocity is as follows. Since travelling bands 
correspond to (staggered) local oscillations, two succes- 
sive passages of bands at a given point of the system 
correspond to the completion of a local oscillation cycle 
for this point. Noting that the condition of a constant (cr) 
imposes that bands disappearing at one end must reap- 
pear immediately at the other, we conclude that each 
band must cross the system in one local oscillation pe- 
riod. Since this is of order rs, and recalling that H = 1, 
we deduce Wband ^ 1/ts for the band velocity. Our data 
show that Wband indeed decreases with rg , and confirm a 
roughly linear trend with t^^ (see Fig. O. 
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FIG. 8: Space-time plots of the stress a{z, t) showing other 
examples of travelling bands, (a) 'Target pattern' analogue 
seen at ts/tm = 30 and (a) = 8.5; (b) Alternating bands 
seen at ts/tm ~ 400 and (a) = 5; (c) Transition between 
oscillating band and travelling band regime seen at ts/tm ~ 
1000 and (a) — 7. Other parameters as previously. 



3. Spatiotemporal rheochaos 

The third main regime encountered in the phase dia- 
gram arises in two disconnected pockets (regions 'C in 
Fig.^ at moderately long values of rg relative to tm and 
for strongly off-centered values of the imposed stress (cr) . 
An increased complexity of oscillations within regime O 
on approaching the wings of the phase diagram was al- 
ready noted previously; in regions C, it finally leads to 
chaotic behavior. 

The spatiotemporal patterns produced in this regime 
are extremely varied and often appear like complex ver- 
sions of periodic patterns that arise in other parts of 
the phase diagram. Examples of rheochaos are given in 
Figure ^1 in (a) , we have a regime of chaotic bands 
with 'random defects' between bands (this is similar to 
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tic phase diagram at lower k might display much larger 
chaotic pockets than Fig. ^ 



VI. LOW-ORDER RESULTS 



FIG. 9: Velocity of travelling bands vs. Tg ^, measured 
in situations where a single band is present in the system: 
(•) {a) = 8.0; (■) (a) = 7.0; (a) (a) = 6.0. The plot shows 
that «band is roughly proportional to Tg^. Error bars are 
±0.15. Parameters as in previous figures. 



rheochaotic patterns observed in shear-thinning micelles 
in Ref. Is^ : in (b), we see a situation which resembles 
the travelling bands of regime T, but here disordered and 
following a somewhat wiggling motion; plot (c) presents 
the case of flip-flopping bands where irregularities ap- 
pear within the bands themselves; plot (d) shows travel- 
ling bands whose nucleation point changes 'at random'; 
in (c) and (f), we observe a 'bubbly' phase of localized, 
short-lived shear bands that appear and disappear errat- 
ically in the cell (although occasionally, a band survives 
and shoots across the system). 

Chaos, for each point so marked in Fig. 0] was con- 
firmed by computing the largest Lyapunov exponent, /iL. 
In Fig. ^2 we give a plot of when (a) is varied on the 
horizontal line ts/tm = 20 of the phase diagram. 

Returning to the data of Fig. ^1 these patterns vary 
in how highly developed is the chaos: for instance, pat- 
tern (b) seems much more erratic than pattern (c) . This 
is reflected on the Lyapunov exponents: we find //l — 1-2 
and /iL — 0.2 respectively. However, we were unable to 
demonstrate true low-dimensional chaos within our high- 
order truncation results. Nor could we numerically ex- 
plore any route to chaos, as the transition window proved 
unattainably small in the phase diagram. We do however 
discuss below the route to chaos within the low-order 
truncation of the model. 

A last comment regards the role played by the value 
of K, the stress diffusivity, in chaotic regimes. Numeri- 
cal solutions in this article were computed with a default 
value of K = 10"^, but data obtained with different val- 
ues suggest that smaller values of the stress diffusivity 
favor chaos. This makes sense, as k is a damping term 
for higher Fourier modes [see Eq. H15|) ]: when k is smaller, 
more modes are involved in the dynamics, and thus the 
system is effectively higher-dimensional. Realistic values 
for K might be much smaller than our choice, which was 
dictated by the need to find A^-independent behavior at 
relatively modest N > 25. Accordingly a more realis- 



The results presented so far all used a high-order trun- 
cation, with a view to obtaining a numerically reliable 
and precise representation of the continuum model of 
Eqs. ©-linjl. In this section, we instead study the be- 
havior in the simplest nontrivial mode decomposition 
of Eqs. 10-®, with A^ = 3. The corresponding low- 
dimensional results provide only a caricature of the full 
dynamics; but they bring interesting insights. This ap- 
proach is similar in spirit to the classical simplification of 
the Rayleigh-Benard equations (describing thermal con- 
vection) into the Lorentz equations [ssf . 

By restricting the number of dynamical variables, it 
permits an analytical or semi-analytical study of some 
of the model's features. It also clarifies which features 
of the model survive in low dimension and thus are in 
some sense robust; we shall see that rheochaos survives 
the truncation, and is not therefore solely the result of 
high-order couplings in the full dynamics. 

Finally, there are interesting physical consequences to 
demonstrating the presence of chaos in a low-dimensional 
Fourier decomposition. The results we have presented so 
far, along with other works p7l Is^. Is^ . might be taken 
to suggest a generic interpretation of rheochaos in com- 
plex fluids as resulting from the erratic motion of discrete 
interfaces between shear bands. This idea is appealing 
because, in turn, it suggests that rheochaos may be more 
efficiently modelled by focussing on interfacial dynamics 
and writing an equation of motion directly for the inter- 
face [23 [53 (rather than for the whole fluid as we are 
doing here). The presence of rheochaos in our low-order 
truncation (which allows only smooth spatial variations 
of a and m) shows that sharply defined bands are not an 
automatic prerequisite for rheochaos. 

As already explained in Sec. IIVI the lowest possible 
mode decomposition allowing chaos is obtained for N = 
3. Since the uniform Fourier modes, ctq = (cr) and mg, are 
dynamically inactive, this reduced representation of the 
model is effectively four-dimensional, with tTi(i), cr2(i), 
TOi (t) and TO2 (t) as the degrees of freedom. Drawing from 
the general analytical expression given in Appendix ^ 
the dynamical equations for the four retained modes have 
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FIG. 10; Various types of spatiotemporal rheochaos observed in the model (high-order truncation). Numerical parameters as 
in Fig. 01 
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FIG. 11: Plot of the largest Lyapunov exponent /il computed 
as a function of the imposed stress (a) on the line ts/tm = 20 
(high-order truncation). Chaotic behavior (positive values) 
appears only on the wings of the diagram. Imprecision on the 
exponent values is typically ±0.1. Numerical parameters as 
in Fig. 01 



the relatively simple form (recall (cr) is a constant): 



CTl 



0-2 



mi 

7712 



-a + 2b{a) - 3c(ay - K{ — 



<J,--ca, 



+ [b ~ 3c((t)] (Ti(T2 — ^ Ami 



^2 



+ - [6 - 3c(fT)] al - -c alo2 - Xm^ 

(cti - mi)/Ts 
(0-2 - m2)/Ts, . 



As with the high-order truncation, we have explored 
numerically the phase diagram for this minimal model. 
We indeed find chaos, confirming that rheochaos is a ro- 
bust feature of the underlying physics. As Figure El 
shows, the global structure of the phase diagram is also 
robust and presents the same features as in the high trun- 
cation: chaotic behavior appears on both wings of the 
diagram (off-centered values of (cr)) and for ts/tm ratios 
not too large (structural relaxation not too slow), while 
periodic states are observed elsewhere. One observation 
is that chaotic pockets are now much smaller than in the 
high-order truncation phase diagram: this is natural as 
chaos is usually facilitated in higher dimensions. 

We do not classify periodic points in this low-order 
model as 'oscillating bands' or 'travelling bands': the 
stress variations are anyway too smooth to allow the def- 
inition of proper shear bands. However, as illustrated in 
Fig. 1131 there remain low-dimensional analogues of the 
oscillating and travelling bands, in roughly the same lo- 
cations as in the high-order truncation. 

In Figure ^1 we show a typical trajectory in the 
chaotic regime of the low-order model: this appears less 
erratic than in the high-order counterpart (Fig. I1(J|) . es- 
sentially amounting to a slightly irregular oscillation of 
two 'bands'. The largest Lyapunov exponent is — 0.4. 



Finally, we studied the route to chaos in the low-order 
truncation, and found a classical period-doubling sce- 
nario (Fig. ll5|l . This result does not, however, imply any- 
thing about the route to chaos in the high-dimensional 
version of the model. 

The results of this section show, not only that sharp in- 
terfaces are not necessary for rheochaos, but also that the 
low-dimensional representation of our shear-thickening 
model can reproduce most of the important physical fea- 
tures observed in higher dimensions. This low-order trun- 
cation could also enable a deeper understanding of the 
oscillating and travelling band regimes: we have not ex- 
plored in detail the relation between rflip, rs, and the 
band velocity. It might also permit a deeper understand- 
ing of the effect of the global stress constraint, partic- 
ularly the mechanism whereby this affect the dynamics 
and makes it more or less complex. 



VII. CONCLUSIONS 

We have introduced a shear-thickening fluid model, al- 
lowing for spatial inhomogeneity, in which the relaxation 
of stress couples to a memory term linked to slow struc- 
tural changes in the fluid. The interplay between the 
rheological and the structural dynamics leads to compli- 
cated spatio-temporal dynamics. The model maps onto 
the FitzHugh-Nagumo model for neural networks, but 
with a different type of global constraint (imposition of 
the average stress (cr)), specific to complex fluids. 

In steady state for both the structure and stress, the 
flow curve is monotonically increasing. However, short- 
term flow curves, valid on timescales too short for the 
structure to relax, are non-monotonic and unstable. This 
obliges the appearance of unsteady shear bands. 

The phase diagram of the model on varying (cr) and the 
ratio of structural and stress relaxation times, ts/tm(> 
4) has a simple overall structure. Three main regimes 
were found. For very large ts/tm, we found oscillating 
shear bands; elsewhere we found periodic travelling shear 
bands at mid-range values of {a) , with regions of spatio- 
temporal rheochaos for off-centre values. 

In the oscillating band regime, shear bands must os- 
cillate synchronously while respecting the imposed value 
for the mean stress (cr). When the volume fractions of 
the two bands within the cell are unequal (off-centre (cr)) 
oscillations become increasingly complex but remain pe- 
riodic. In the other periodic regime seen in the model, 
bands nucleate at a border point of the container, or 
more rarely at an interior point, and travel across it at a 
near-constant velocity which scales roughly as . 

Our model exhibits spatio-temporal chaos. The corre- 
sponding space-time patterns are varied, but fall into a 
common picture (coherent with that of other groups js^ 
Is^): rheochaos manifests as a flow that restlessly at- 
tempts to form steady shear bands, but fails due to in- 
ternal structural constraints. 

It is known that the constraint on (cr) (set by the 
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FIG. 12: Phase diagram in the low-dimensional truncation of the model: (•) Periodic points; (a) Chaotic points. Note that, for 
pictural clarity, a horizontal spacing between points A(a) = 0.25 was chosen for the plot, but our actual numerical exploration 
of the phase diagram was performed on a much finer grid with A{a) = 10~^. The outer dashed line delineates the limit of 
linear stability [Eq. 0]. Dotted lines around chaotic pockets are merely guides to the eye. Numerical parameters as in Fig. El 
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FIG. 13: Low-dimensional analogues to the oscillating band 
regime (a) and to the travelling band regime (b) . Note the ill- 
defined 'interface' between 'bands'. Parameters: (a) ts/tm = 
Iff*, (a) — 7; (b) ts/tm ~ 40, (a) — 7; other parameters as in 
Fig.H 
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FIG. 14: Rheochaos in the low-dimensional representation 
of the model, here plotted for ts/tm = 60 and (a) — 3.55. 
Numerical parameters as in Fig. |^ 

mechanical torque on a Couette) is essential for shear- 
thickening fluids as it sets the respective fractions of the 
high- and low-shear bands in the flow '^l ; and that flnite 
stress diffusivity k provides a selection criterion among 
possible banded flows 0. Within our model, the im- 
portance of these quantities extends to time-dependent 
cases: dynamical complexity increases as (a) moves off- 
centre, whereas a smaller diffusivity k promotes chaos. 

Taking advantage of our Fourier-space representation, 
we also studied a four-dimensional truncation for which 
we showed that the important physical features found in 
high dimension persist, including the general structure of 
the phase diagram. 

The presence of chaos in the low-order model shows 
that sharp interfaces between bands are not necessary 
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FIG. 15: Route to chaos in the low-dimensional truncation, 
shown for ts/tm = 60. Plotted are three-dimensional pro- 
jection on the ((Ji, a2, mi)-space of the four-dimensional at- 
tractors. From (a) to (d), first steps of the period-doubling 
cascade: (a) (a) — 3.53, (b) (a) = 3.535, (c) (a) — 3.5375, 
(a) (a) = 3.5379; (e) Chaotic attractor obtained for (a) = 
3.55, corresponding to the situation shown in Fig. 1141 Nu- 
merical parameters as in Fig.|3] 

for rheochaos. This does not mean that such interfaces 
have no effect, only that the chaos would persist without 
them. 

In conclusion then, our work, along with others' [s^ . 
Is^ . supports a 'frustrated shear-banding' picture of 
rheochaos, in which flow inhomogeneity (whether sharp 
interfaces or smooth variations) plays a crucial role. A 
gap between theory and experiment still remains: our 
model makes no direct connection with the microscopies 
of the considered fluids. Improved experimental infor- 
mation about the dimensionality of the observed chaos, 
or the route into, for shear-thickening systems would be 
most helpful. One general prediction of our model that 



may be easily tested experimentally is that the dynam- 
ics should get more complex as the imposed stress (tr) is 
chosen closer to the edges of the fluid's unstable window. 

One potentially important factor that has been ne- 
glected in the present model is the role played by normal 
stresses in the fluid. In many of the experimental fluids 
we are concerned with, their magnitude is significant in 
high-shear regimes. One could also consider alternative 
scenarios for rheochaos, where stress variations do not 
couple to nonconserved structural changes as proposed 
here (these include modulation of concentration at finite 
wavevector, i.e. microphase separation) but instead to 
conserved modes such as global concentration fields. 
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APPENDIX A: NON-DIMENSIONALISATION 

We here describe how the non-dimensional equa- 
tions for the model are obtained from the original 
physical equations. This is especially useful if one wants 
to deduce the numerical value of a non-dimensional pa- 
rameter used in the model from an experimental value. 
Throughout this Appendix, we use the convention that 
normal letters are for dimensional, physical quantities, 
and starred letters for non-dimensional ones. In the rest 
of the article, all quantities are non-dimensional by de- 
fault, so that stars have systematically been omitted. 

The dimensional version of Eq. writes 

a = Gj- Tj^i V + - ca^ (Al) 
-A J Tg-ig^*-*')/"^ ait') dt' + kV> 

where G is the transient elastic modulus, and other quan- 
tities have been defined in the main text. 

We use G, tq (defined below) and H (the Couette axial 
extent) as, respectively, the units for stress, time and 
space, and define the corresponding reduced quantities: 
a* = a/G, t* = t/ro and z* = z/H. 

Substituting these quantities in Eq. IjAip . we obtain 
the non-dimensional expression 

cr* = 7 _ 4- 6G'Toa*2 - cGVocr*^ (A2) 

We now have to choose a value for the unit time tq: 
for practicality, because we always work on timescales 
longer than tm, we choose tq = IOOtm. We next define 
the reduced Maxwell time = tm/tq and the constant 
a*, its inverse: we have rj^ = 0.01 and a* = t^"^ = 100. 
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Finally, we transform the remaining physical param- 
eters of the model into reduced quantities: b* = ^Ctq, 
c* = cG^tq, X* — Xtq, Tg — ts/tq, k* — ktq/H^. In 
terms of these reduced quantities, Equation IA2I rewrites 

a* = 7 - a*o-* + 6*0-*^ - c*(T*^ 

-\*J T*-ie(*'-*'*)/^s a*{t'*) dt'* + K*V^,a* 

which is exactly the expression given in Eqs. CQ-Q of 
the main text (dropping stars out). 

APPENDIX B: FULL ANALYTICAL 
EXPRESSION FOR MODE EQUATIONS 

In this Appendix, we give the full analytical expression 
of the evolution equations for the stress modes cr„(t), of 
which only a shortened version was given in Eq. H15|l . [For 
memory modes, the analytical expression is straightfor- 
ward and was given in Eq. I|16(l .] 

I 



-(3- 



where S is the usual Kronecker symbol, and the following 
shorthands were used 

qi = max{n + p — N + 1,0} 

q2 = N -l + inm{n+p- N + 1,0} 

qs = max{0,j3— n} 

54 — — 1 + min{0,p — n}. 

In high-order truncations (N ~ 40), a typical mode 
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